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The External Verification Analysis approach to code verification is extended to solve 
the three-dimensional Navier-Stokes equations with constant properties, and is used to 
verify a high-order computational aeroacoustics (CAA) code. After a brief review of the 
relevant literature, the details of the EVA approach are presented and compared to the 
similar Method of Manufactured Solutions (MMS). Pseudocode representations of EVA’s 
algorithms are included, along with the recurrence relations needed to construct the EVA 
solution. The code verification results show that EVA was able to convincingly verify a 
high-order, viscous CAA code without the addition of MMS-style source terms, or any 
other modifications to the code. 


Nomenclature 


x, y , 2 Cartesian coordinates 

t Time 

u, u Exact and manufactured solution to the linear advection equation 

/, g Arbitrary differentiable functions 

a Advection speed for the linear advection equation 

b Arbitrary constant used in u 

S Linear advection manufactured solution source term 

p Density 

a 1/p 

u x-velocity 

v y-velocity 

w ^-velocity 

p Pressure 

<t>, 0 Navier-Stokes mechanical and conductive viscous terms, respectively 

X x, y, or z 

U a, u, v, w , or p 

7 Ratio of specific heats 

R Ideal gas constant 

k Thermal conductivity 

p, Dynamic viscosity 

e Error norm 

d>, 0 Reference and PDE code solutions, respectively 

h Discretization measure 

A, B h-independent coefficients 

p Observed order-of-accuracy 

c t> Primitive flow variable initial condition (a, u, v, w , p) 
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Xo, 2/o, Zo 
fox , ky , kz , 
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(kAx)* 
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Constants in 4> 

Constants in <j> 
constants in <f> 

Wavenumber of a single Fourier component 
Numerical wavenumber of a spatial differencing scheme 
Spatial differencing scheme’s relative error 
Frequency of a single Fourier component 
Number of grid points in test case grid 
Time marching scheme amplification factor 
Local and global time-marching error, respectively 


I. Introduction 


As computational aeroacoustics codes become more sophisticated and widespread, the need for quanti- 
fying their accuracy becomes more important. The field of Verification & Validation (V&V, cf. the recent 
book 1 on the subject) is focused on developing techniques for assessing the correctness and reliability of 
scientific codes, including determining the accuracy of their numerical algorithms (i.e. , verification) , and the 
closeness with which their output correspond to nature or experiment (i.e., validation). 

Code verification, the crucial first step of V&V, is an investigation into whether a CAA code (or any 
PDE solver) is “solving the equations right.®’ Specifically, a code can be deemed verified if it is rigorously 
shown to solve its intended PDE to the formal order-of-accuracy of its numerical schemes. Code verification 
is distinct from solution verification, which is the process of estimating the amount of numerical error in a 
PDE code’s solution. Se£E and especially 4 for clear discussions of code verification, solution verification, 
and other V&V concepts. 

Code verification requires evaluating the error in a PDE code’s solution, which in turn requires the 
comparison of a PDE code’s output to a reference solution. One approach to obtaining such a reference 
solution is called the method of exact solutions (MES, a term likely coined by Knupp and SalarP). One 
simply finds an exact solution to the PDE the code in question solves somewhere in the literature, applies 
the PDE code to the problem to which the exact corresponds, and compares the code’s output to the exact 
solution. For example, to verify a PDE code that solves the linear advection equation, 


du du 

m +a d^ = 0 ' 


(i) 


the exact solution 

u(x,t) = f(x - at), (2) 

where / is any differentiable function, is taken as the reference solution and used to calculate the PDE code’s 
error. 

Equation ([2]) is a quite general solution to ([!]): one has considerable flexibility in choosing /, and thus the 
form of the solution. Unfortunately, known exact solutions to non-trivial PDEs are usually not as flexible, 
instead tending to be either too simple to adequately exercise all the capabilities of the code, or are tedious to 
evaluate numerically. An example of an exact solution that is simultaneously too simple and too complex is 
provided by Knupp and Salari,®who reproduce the exact solution to unsteady homogeneous heat conduction 
in a sphere. Evaluating the solution involves two infinite series, a triple integral, and various flavors of Bessel 
functions, Bessel function roots, and Legendre polynomials. Such a solution, however, would be inadequate 
for verifying a general heat conduction code that can handle arbitrary geometries and non-constant physical 
properties, since it (the solution) does not exhibit any of these characteristics. 

The Method of Manufactured Solutions (MMS) was developed to overcome this limitation of MES. First 
introduced by Steinburg and Roache,® MMS can be used to make nearly any expression (within reason) an 
exact solution to a code’s PDE at the cost of introducing an additional source term to PDE, and thus the 
code itself. Excellent explanations of the MMS technique are available in the hteratureJ390® but a brief 
example of its use is given here. Suppose, again, that a PDE code solves the linear advection equation 0 . 
and one wishes to use MMS to verify the code. The reference ( “manufactured” ) solution is specified next — 
or, as Roache® points out, before even looking at the PDE. A simple choice is made 


u{x, t ) = sin(x — bt) 


( 3 ) 



American Institute of Aeronautics and Astronautics 


where u is the reference solution and b an arbitrary constant. The proposed solution is then substituted into 
the governing equation to find the MMS source term: 


S(x,t) = [sin(x — bt)] t + a [sin(x — bt)\ x 
= —b cos(x — bt) + a cos(x — bt) 

= (a — b) cos{x — bt) . (4) 


The source term S(x,t) is finally added to the original governing equation: 

du du . ... 

— — h a— = b( x, t) = (a — b) cos(x — bt). (5) 

dt ox 

The manufactured solution u = sin(x — bt) is an exact solution of the modified equation ([5]), and can be used 
as a reference solution when verifying any PDE code that can be made to solve ([5|. (Of course, if one chose 
b = a, then S(x,t) = 0 and u is an exact solution to the original PDE. Doing the same for more complex 
PDEs is generally not possible.) 

The Method of Manufactured Solutions has enjoyed fairly wide adoption in the CFD community by those 
concerned with V&V, having been used to verify Finite Volume method CFD codes!®®!] solving both the 
Euler and Navier-Stokes equations on unstructured meshes, fluid-structure interaction (FSI) codesP^EED anc [ 
RANS codes with turbulence models. 17 !® Instances of its use in non-CFD applications are less prevalent, but 
do exist. Merroun et al. 19 showed how MMS could be used with a code simulating convective heat transfer 
in multiple channels (the code’s application was nuclear reactor cooling), and Pautz 20 communicated results 
of applying MMS to a finite element code solving equations with applications to nuclear transport. Also, 
one of the FSI examples previously mentioned® had its origin in biomechanics problems. 

Despite its generality and many proponents in the V&V community, MMS appears to be used to verify 
new and existing PDE codes infrequently. For example, a review of recent issues of the Journal of Com- 
putational Physics (JCP, issue 274) and Computers & Fluids (C+F, issue 103) was performed, and the 
code verification technique(s) used by the authors, if any, were noted. Results are summarized in Figure [l] 
where the number of papers using MES and MMS are displayed, along with Code Comparison,^ a less- 
rigorous code verification technique that involves comparing the output of one code with another. Clearly, 
the Method of Manifactured Solutions is the least-popular code verification technique. The reason for this 
is likely related to the MMS source term that must be added to the PDE code, which does tend to be com- 
plex for non-trivial governing equations (though computer-algebra systems like Mathematica, Matlab make 
this complexity much more manageable). The ideal code verification tool would allow the user to create a 
reference solution consisting of elementary, easy-to-evaluate functions, but would not require any changes 
to the equations the PDE code solves — the flexibility of MMS with the unobtrusiveness of MES. External 
Verification AnalysiP^l i s the authors’ attempt at such a tool. 


II. External Verification Analysis 


External Verification Analysis’ (EVA) point of departure is related to, but different from, the Method 
of Manufactured Solutions. Instead of specifying the complete solution to a PDE up front (i.e., the MMS 
approach), a solution on a reduced-dimensional surface is chosen, and then expanded beyond the initial 
surface using a Taylor series. The coefficients of the Taylor series are calculated using the PDE itself and 
Cauchy-Kowalewski recursion®®! For example, with the same PDE from the MMS demonstration, the 
solution on the t = 0 “surface” is specified, 


u( x,0) = u| 0 (x). 


( 6 ) 


If a solution similar to that of the MMS demonstration is desired, u| 0 (x) = sin(x). The solution for t ^ 0 
is found using the Taylor series 


u(x, t) 


u\ 0 +t 



t 2 d 2 u t 3 d 3 u 
2 dt 2 0 6 dt 3 


( 7 ) 


Expressions for the temporal derivatives in the coefficients of 0 are found by performing Cauchy- 
Kowalewski (CK) recursion on the linear advection PDE. First, the temporal derivative is isolated: 


du du 

dt a dx 


( 8 ) 
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Figure 1. Counts of code verification techniques used in recent issues of the Journal of Computational Physics 
(JCP) and Computers & Fluids (C+F). 


which, after substituting t = 0, provides an expression for the first temporal derivative in the Taylor series, 
i.e., 


du 

~di 


du 

dx 


( 9 ) 


where jt^ | Q is found by differentiating ©■ The second-order temporal derivative is found by first differenti- 
ating Q with respect to t: 

d 2 u d 2 u 

dt 2 = ~ a dx dt 


An expression for the mixed spatial-temporal derivative g x g t 
differentiating with respect to x , i.e.: 

d 2 u d 2 u 

= — a 


in ( 10 ) is found again from 


( 10 ) 

this time by 


( 11 ) 


dx dt dx 2 ' 

— i — a 2 a 2 

Taken together, (10 1 and (11) provide the necessary information to calculate (Ay at t = 0: first is found 


by differentiating ([6]) twice, which is used to evaluate ® x g t from (IT] ), which in turn allows one to find 
from ( 10 ). This recursive process of successively finding expressions for temporal and mixed spatial-temporal 
derivatives in terms of pure spatial derivatives can be repeated indefinitely, allowing the Taylor series of ([T]) 
to be extended to an arbitrary order. 

The preceding discussion shows that there are essentially two phases to CK recursion: first, an expression 
for a needed derivative of u is found by differentiating the governing equation; second, the temporal order of 
the terms in the new derivative expression are reduced, again by differentiating the governing equation. This 
process is described in pseudocode as the recursive procedure CK in Algorithm [lj The procedure takes two 
arguments: derivdd, some object that represents the order and dependent variable that should be found (i.e., 
the programming equivalent of a symbol like fK T, ) , and deriv .cache, an object that stores the derivatives 
that have already been found and their numerical values. The routine begins by checking if the value of the 
derivative derivdd has been calculated previously and is thus stored inside deriv.cache. If it has, the routine 
does nothing (pass), the rest of the if block is skipped, and the program exits the routine. 
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Algorithm 1 Pseudocode for Cauchy-Kowalewski recursion 
procedure CK(deriv_id, deriv_cache) 
if deriv_id in deriv .cache then 
pass 

else if GET_TORDER(deriv_id) == 0 then 
deriv _val t— EVAL_lc(deriv_id, deriv .cache) 

STORE_VAL(deriv_id, deriv.val, deriv .cache) 
else 

expr 4 — GET_EXPR(deriv_id) 
for all local deriv id in expr do 

CK(local_derivJd, deriv.cache) > Recursion! 

end for 

deriv _val t— EVAL_EXPR(expr, deriv.cache) 

STORE_VAL(deriv_id, deriv.val, deriv.cache) 

end if 

end procedure 


If derivdd does need to be found, the CK routine determines if the needed derivative is a pure-spatial 
derivative (e.g., |A() or a mixed spatial-temporal or pure-temporal derivative (e.g., -§^rgi or §jr), repre- 
sented in the pseudocode as the check in the GET_TORDER(deriv_id) == 0 part of the else if statement. 
GET.TORDER is imagined to be a function that takes the deriv _id object and returns the order of its temporal 
derivative (so GET_torder( g c \ u at ) is 1 and GET_torder( ) is 3). If derivdd is indeed a pure-spatial 
derivative, its numerical value is calculated by the EVAL.IC function and then saved in the deriv.cache object 
by the STORE.VAR routine, and the routine returns. 

The last branch of the if / else if / else block handles the case when derivdd is a mixed spatial-temporal 
or pure-temporal derivative than has not been calculated previously. First, an expression for the desired 
derivative is found by differentiating the governing equation and placed in the expr variable, represented in 
Algorithm l] by the line containing the GET.EXPR function, which is the pseudocode equivalent of obtaining 
(10 1 or (Jll) through the differentiation of ([8]). Next, each term in expr is passed to the CK routine, “winding 
up” the recursion. The recursive CK call has the effect of finding any unknown derivatives in expr and storing 
them in deriv.cache. Finally, the “unwinding” of the CK procedure is accomplished by the last two lines of 
the else branch, which find the numerical value of the derivdd derivative using the previously-calculated 
derivative data in deriv.cache and the EVAL.EXPR routine, and then save the value in the deriv.cache variable, 
again using the STORE _VAL routine. 

Algorithm [l] could be used to easily implement Cauchy-Kowalewski recursion as part of an EVA tool, 
especially if a computer algebra system (CAS) is available. Indeed, such a program has been devolved using 
the Python programming language 25 and SymPyJ^ an open-source Python CAS library, and subsequently 
used to very the EVA tool itself. While the flexibility provided by a high-level language like Python and 
CAS library like SymPy is attractive, such flexibility comes at the cost of efficiency. Also, unless a compiler 
optimizes it, recursive routines tend to be considerably slower than an explicit loop-based representation. 
For these reasons, unwinding the recursion in Algorithm [l] and finding an expression for the arbitrary-order 
derivative of the targeted governing equations will be discussed in the next section. 


III. EVA with the Navier-Stokes Equations 

This work proposes to use the technique presented in the previous section to construct a reference solution 
to the Navier-Stokes equations with constant properties (i.e., ratio of specific heats, gas constant, viscosity, 
and conductivity) that can then be used to verify a high-order viscous CAA code. The three-dimensional 
form of the Navier-Stokes equations with primitive variables a = 1/ p, u, v, w , and p are used with the EVA 
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tool: 
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(12a) 

(12b) 

(12c) 

(12d) 

(12e) 

(13) 


$ = p 


du du 
dx dx 


dv dv dw dw du dv du dw dv dw 
dy dy dz dz dx dy dx dz dy dz _ 

dv dv dw dw du du dw dw du du dv dv 
+ dx dx + dx dx ^ dy dy + dy dy dz dz + dz dz 


+2 


dv du dw dv du dw 
dx dy dy dz + dz dx 


(14) 


Again, ( 12 )-( 14 ) assume the dynamic viscosity p , thermal conductivity k, and gas constant R are all constant. 

As discussed in the previous section, an explicit expression for the arbitrary-order derivative of the 
governing equations will allow for a more efficient implementation of the EVA technique. Because each 


equation in ( 12 ) only contains sums of products, all that is needed is nested applications of the Leibniz rule 

d n [f{x) ■ s(z)] 


dx r 


= £ 

ra— 0 


n 


l f d m g 


mj dx n ~ m dx m 


(15) 
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Applying (15) to the first equation in ( 12a| ) (the continuity equation) gives 

da da 


da ( da 
= ~ I u — 
dx 
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dx b ~ a dy d ~ c dzf~ e dt n ~ m dx a+1 dy c dz e dt m 

f^b—a-\-d—c-\-f—e-\-n—my ^a+c+l+e+m^- 

+ dx b ~ a dy d ~ c dzf~ e dt n - m dx a dy c+1 dz e dt m 

Qb—a-\-d—c-\-f—e-\-n—myj ^a+c+e+l+m^. 


(16a) 


(16b) 


n 
x ra 

m = 0 
^a+l+c+e+m^ 


dx b ~ a dy d ~ c dzf~ e dt n ~ 1 

^b—a-\-d—c-\-f—e-\-n—m^j 

dx b ~ a dy d ~ c dzf~ e dt n ~ 1 

ga+l+c+e+m u 

dx a+1 dy c dz e dt m 

ga+c+l+e+m 

dx a dy c+1 dz e dt m 

ga+c+e+l+m w 

dx a dy c dz e+1 dt n 


dx a dy c dz e+1 dt r ‘ 


(16c) 


The continuity equation is rearranged to isolate the temporal term in (16a). Equation (16b) is an ex- 


pression for the ?rth temporal derivative of the continuity equation, found by applying the Leibniz rule to 


f-differentiation and (16a). Finally, the Leibniz rule is again used (three times) with (16b) to obtain (16c), 
the full expression for any x, y , z , t derivative of the continuity equation. The corresponding expressions 


for the other equations in (12) are given in the Appendix as (50) (55). Taken together, (16c) and (50) (55) 
are an explicit representation of the GET.EXPR and EVAL_EXPR functions in Algorithm [1] i.e., they allow 
EVAL_EXPR(GET_EXPR(deriv_id), deriv_cache) to be implemented in a computer code without resorting to a 
CAS. 


Equations ( 16c ) and ( 50 1 ( 55 ) are called recurrence relations because they express higher-order temporal 


derivatives of the Navier-Stokes equations in terms of lower-order derivatives. Understanding how this 


recursive process works is made easier by expressing ( 12 ) in the more-general form 

dU f dU dU dU d 2 U d 2 U d 2 U d 2 U d 2 U d 2 U 
~dt ~ i V’ dx ’ dy ’ ~dz ’ ~d^ 1 dxdy' dx dz ’ ~d^ ’ dydz’ dz 2 ’ 


(17) 


where U are the unknowns <r, u, v, w, p. The spatial coordinates are collapsed into one unknown in the 
expression 

dU 


dt 


= f[U, 


dU d 2 U 
dX’dX 2 ’ 


(18) 


where the notation stands for any derivative with respect to x a y b z c for all a + b + c = n (so, for 

example, indicates the derivatives §^, and g^i). Next, consider each term 

on the right-hand side of (18). When building the Taylor series in ,0, the derivatives will naturally be 
calculated in increasing order, i.e., first U will be found, then ggg etc.. If this is so, then a value for 
U will already be known when attention is turned to calculating but gg an d §yt will not. To express 
this dependency, the notation 

dU dU d 2 U 


dt 


dX ’ dX 2 
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will be used, and means and §y? are needed to calculate Equation (19) is expressed graphically 
in Figure [2j 



Figure 2. 


CK recursion dependency diagram for ( |19| ). 


Next, consider the equivalent of ( 18 ) for 


d 2 u 

dt 2 


After recognizing that all terms in ( 12 )-( 14 ) are sums of 


products, one sees that the correct relationship is 


d 2 U / dU dU d 2 U d 2 U d 3 U \ 

~W v ’ ~dt' ax 7 ’ dxdt’dx*’ dx 2 at') ' 


(20) 


Compare (20) to ( [l8| ). If only the new terms in the right-hand side of (20) are retained (i.e. , those that are 
found in (20) but not (|l8|), then the resulting dependency expression for is 


d 2 U d 2 U d 3 U 
~W dXdt' dX 2 dt' 


( 21 ) 


which is recognized as the same expression that would have been found by differentiating each term in 
the right-hand side of ( 19 ) with respect to t. Knowing this, a general dependency relationship for the 
temporal and mixed spatial-temporal derivatives of the Navier-Stokes equations can quickly be determined 
by differentiating (19) with respect to t n and X a : 


Qa+n+ljj ga+l+njj Qa+2+njj 

dX a dt n+1 dX a+1 dt n ’ dX a + 2 dt n 


(22) 


for a, n > 0. 

Equation (22) allows one 
tives. For example, Figure ]3 
derivatives are left; Figure 4 


to create graphical dependency diagrams like Figure [2] for higher-order deriva- 
shows the result of beginning with and using ( |22| until only pure-spatial 
shows the same for . 



Figures [2]-[4] are top-down representations of the Cauchy-Kowalewski process — they begin with 
desired result (^r) and end with the input, or starting point { g X m )• One would prefer to unwind 


the 

the 
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• • • q3 t t ( 

Figure 4. CK recursion dependency diagram for and its dependents. 


recursion and instead start from the new pure spatial derivatives needed for a given temporal derivative, 
then find the required mixed spatial-temporal derivatives, then finally the desired pure-temporal derivative. 
This can be achieved by simply eliminating all but the right two nodes of each row in Figures [2] |4j and then 
moving from left-to-right, bottom-to-top with those that remain. After removing the unnecessary nodes, 
Figure [3] becomes Figure |5j and Figure [2J Figure [6] (Figure [ 2 ] is unchanged). 



The “left-to-right, bottom-to-top” calculation procedure just described is represented in pseudocode in 
Algorithm [2] Here, the deriv_cache variable is imagined to be some object that holds the numerical values 
of derivatives calculated previously (identical to its meaning in Algorithm [lj, and tot_order is an integer > 0 
equal to the order of the pure-temporal derivative that will be calculated by the CK.UNWOUND procedure 
(e.g., tot_order = 2 for Figure [5] and tot_order = 3 for Figure [h]). X.order and t.order are also integer 
variables, indicating the spatial and temporal order, respectively, of the derivatives that will be calculated 
by the EVAL_DERiv routine and stored in the deriv .cache variable — for example, if X.order = 2 and 
t.order = 1 are passed to EVAL.DERIV, then g d x l > g t f° un d and saved in deriv.caclie. 

A pseudocode listing of the EVAL.DERIV routine is shown in Algorithm [3j The outer if /else block 

9 ofl34l 
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Figure 6. CK recursion dependency diagram for - and its dependents with unnecessary nodes removed. 


Algorithm 2 Pseudocode for unwound Cauchy-Kowalewski recursion with the Navier-Stokes equations, 
procedure CK_UNWOUND(tot_order, deriv .cache) 

X.order 4 — 2 • tot.order 
for t.order 4 — 0, tot.order — 1 do 
X.order 4 — X.order — 1 
EVAL_DERlv(X_order, t.order, deriv.cache) 

X order 4 — X order + 1 
EVAL_DERlv(X_order, t.order, deriv.cache) 

X .order 4 — X.order — 2 
end for 

EVAL_DERlv(X_order, t order, deriv cache) > Here X.order = 0, t order = tot .order 

end procedure 
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decides if the desired derivatives are temporal/mixed-spatial-temporal or pure spatial derivatives of the initial 
condition. The code for the two branches is similar: both contain loops that iterate over the derivatives 
indicated by the X.order and t.order derivatives (and, therefore, the g Xa notation), and both eventually 
use the STORE.VAL routine to save each iteration’s work in the deriv.cache variable. The x.order, y_order, 
and z_order variables are integers that correspond to the x-, y. and z-order of the derivative calculated by 
either EVAL_RECUR_REL (using the Navier-Stokes recurrence relations in (16c) and (50) (55)) or EVAL.IC (by 
differentiating the known initial condition). 


Algorithm 3 Pseudocode for EVAL.DERIV in Algorithm [2j 
procedure EVAL_DERiv(X_order, t_order, deriv_cache) 
if t order > 0 then 

for z order <— 0, X order do 

for y .order 0,X_order — z.order do 
x.order X.order — z.order — y.order 

deriv.val EVAL_RECUR_REL(x_order, y.order, z.order, t.order, deriv .cache) 
STORE_VAL(x_order, y.order, z.order, t.order, deriv.val, deriv .cache) 

end for 
end for 
else 

for z.order 0, X.order do 

for y.order 0, X.order — z.order do 
x.order X.order — z.order — y.order 
deriv.val 4 — EVAL.IC (x.order, y.order, z.order, deriv.cache) 
STORE_VAL(x_order, y.order, z.order, t.order, deriv.val, deriv.cache) 
end for 
end for 
end if 

end procedure 


Taken together, the recurrence relations (16c), ( 50 )— ( 55 ) and the pseudocode in Algorithms [2] and [3] can 
be used to efficiently implement Cauchy-Kowalewski recursion and build a Taylor series solution 0 to the 
Navier-Stokes equations. The Taylor series will not, of course, contain an infinite number of terms, and 
thus involve some truncation error. In this work, an estimate of the truncation error is obtained by adding 
together the magnitude of the last two Taylor series terms, and then compared to a desired error tolerance 
set by the user of the EVA tool. If the error estimate is less than the user-specified value, the order of the 
Taylor series is increased; otherwise, the EVA solution for the current spatial location is deemed adequate 
and attention is turned to the next grid point. 


IV. Results 


A. BASS 

The CAA code verified in this work is NASA Glenn’s Broadband Aeroacoustic Stator Simulator 27 (BASS), a 
high-order, block-structured finite-difference LES solver. BASS is capable of working with the fully nonlinear, 
Euler or Navier-Stokes equations in conservative form and two or three dimensions — three-dimensional 
Navier-Stokes results will be presented here. 

A variety of time-marching and spatial differencing schemes are implemented in BASS. The spatial stencils 
available include the Dispersion Relation Preserving (DRP) scheme of Tam and WebfpS with coefficients 
from Tam and Shenj25 and Hixon’s prefactored form 30 of the compact 6 th -order scheme of Lele 3 ^ (C.6), 
as well as standard 2 nd - and 6 th -order explicit central differencing (E.2 and E.6, respectively). For time- 
marching, the four- and five-stage Runge-Kutta schemes of Jameson^ (RK4L and RK5L, respectively), the 
5/6-stage two-step Runge-Kutta scheme of Stanescu and Habashi^ (RK56), and the High- Accuracy Large- 
step Explicit Runge-Kutta (HALE-RK) 7-stage and 6/7-stage two-step schemes of Allampalli et alP^ (RK7S 
and RK67) are used in this work. 
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B. Verification Procedure 

The procedure used to verify the BASS code with EVA is summarized in the following steps: 

1. Set up the problem - specify the initial condition, final time level, extent of the computational domain, 
etc., 

2. Use EVA to find the solution at the final time level, 

3. Use BASS to find the solution at the final time level with a range of discretization measures (grid 
spacings for spatial verification, time step sizes for temporal), 

4. Calculate the error in the BASS solutions through comparison with the EVA solution, 

5. Calculate the error convergence rate (or observed order-of-accuracy). 

In reference to step [4j two different error norms were in this work: the / 2 


e h 



(23) 


and the /max 


ez 




= max 

n 



(24) 


The notation *I> and </>!*) indicate the EVA and i-th BASS code solutions, respectively, with the subscripts 
representing the solution at each of the N spatial locations in common between <f> and (j>. The Z max norm 
tends to be more volatile than the / 2 - 

In step [5] the error convergence rate was calculated using two different methods. The first approach 
assumed the BASS solution error could be well-approximated by the expression 


e* = <t>i - $ 


AhP 


(25) 


To find the convergence rate p, the ratio of the error of two PDE code solutions is formed 


U-i 


Ah? 

AhU 


h? 

hU 


1 


(26) 


and then rearranged to isolate p, 


log (e,) - log (ej_ 1 ) 


lQ g(uu) 

l og ^_Ai_^ log (hi) - log (hi-!) ' 


(27) 


The second expression for p in (27) shows the connection between the convergence rate of the error and the 
error itself: the former is the slope of the latter on a log- log plot. 

The second convergence rate calculation method makes a slightly more general assumption about the 
BASS error, specifically 

ei » Ah? + B, (28) 

where B is a part of the error insensitive to the discretization measure being refined (grid spacing or time 
step). Equation 28 contains three unknowns (A, B : and p), and each e-h combination is associated with one 
BASS code calculation — thus error norm data from three runs will be needed to solve for the unknowns in 
©, i.e., 

Co = A/iq + B 


d = Ahf + B 
£2 — Ah^ + B 1 


( 29 ) 
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which is a non-linear system with no explicit solution. Following Oberkampf and Roy,® the unknowns A and 
B may be eliminated with a bit of manipulation, yielding 


K - h\ _ e 0 - ei 
K ~ K ei - e 2 ’ 


(30) 


one equation for the unknown p. Unless (as in Oberkampf and Roy) hi = r l ho , Equation 
solved explicitly and a root-finding algorithm must be used to find p — the approach taken 


(30) cannot be 
here. 


C. Test Case Parameters 

For the present test cases, the initial distribution of the primitive flow variables were of the form 


(t>{x,y,z) =<p + </>exp ( - 
sin ( k x x 


In 2 


0.15 2 L 
k y y + k z z 


(x - xq ) 2 + {y- y 0 f + {z- z 0 f 

0 ), 


(31) 


with the values of the various constants given in Table [l] The location of the Gaussian perturbations peaks 
(i.e., Xq, 2/0; ~o from ( |31| ) ) for each of the flow variables were placed at different locations in the domain to 
avoid any symmetry in the initial flow; Figure [ 7 ] uses colored spheres to represent these location. The values 
of k x , k y , and k z were set such that the direction the sin() term varies the most rapidly differs for each 
variable, but has the same wavelength of 0.25, giving an approximate wavenumber k = = 8n. Figures^ 

and [9] show slices of the density and ^-momentum initial conditions, respectively. 


0 

(t> 

0 

x 0 

2/o 

zo 

fox 

ky 

k z 

0 

a 

1. 

0.001 

-0.1 

-0.1 

-0.1 

0.875917 

0.0459049 

25.1174 

7.5 

u 

0.03 

0.006 

-0.05 

0.05 

0.05 

4.00477 

14.9460 

19.8048 

9. 

V 

0.02 

0.004 

0.05 

-0.05 

-0.05 

-20.2615 

13.1580 

6.92752 

10.5 

w 

0.01 

0.005 

-0.05 

0.05 

-0.05 

-18.3538 

-14.8626 

-8.59590 

12. 

p 

1 

7 

0.01 

0.1 

0.1 

0.1 

5.03652 

-13.1205 

-20.8359 

13.5 


Table 1. Values of constants for the initial condition found in •mi- 


The common ratio of specific heats value for air, 7 = 1.4, was used, which, when combined with the 
choices p = ^ and p = 1 give a far-held acoustic velocity of unity. The maximum propagation speed 

based on far-held quantities, then, would be c = 1 + \/0.03 2 + 0.02 2 + 0.01 2 ~ 1.0374. A dynamic viscosity 
/i = 0.0002 was chosen, which, when combined with a reference length L = 2 • 0.15 (twice the Gaussian half- 
width in the initial condition), V = -pO.03 2 + 0.02 2 + 0.01 2 (far-held velocity), and p = 1, gives a Reynolds 
number _ 

Re = « 56.1, (32) 

fl 

much lower than what is typically found in aerospace applications. The rationale behind using such a small 
value is related to the usual physical interpretation of the Reynolds number: since it represents the ratio of 
the inviscid and viscous terms, a Reynolds number closer to unity indicates these terms are approximately 
the same magnitude, and problems with one is unlikely to be hidden by the other. Likewise, the Prandtl 
number Pr 1 often described as the ratio of the viscous and conductive terms, was set to 0.7 (the typical 
value for air). Finally, the Peclet number, dehned Pe = Re ■ Pr and a measure of the ratio of the inviscid 
and conductive terms, was Pe = 56.1 • 0.7 = 39.3. 

Because EVA currently solves the Navier-Stokes equations on an inhnite domain, the enforcement of a 
particular boundary condition along surfaces in the solution domain is not possible. BASS, like all PDE 
codes, requires boundary conditions on all domain boundaries, however, and this represents a potential 
discrepancy between the EVA and BASS solutions. To minimize the possibility of mismatched flow solutions 
at the BASS domain boundaries from negatively impacting the verification results, the Giles non-reflecting 
boundary conditiorPS and an initial condition with a limited spatial extent (see Figure [7]) was used. 
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Figure 7. Outline of the BASS verification test case grid with spheres marking the peaks of the Gaussian 
perturbations of the initial flow, with purple, red, blue, orange, and green indicating the cr, u , v, w , and p 
centers, respectively. The radius of each sphere is 0.15, the half-width of the Gaussians. 



Figure 8. Slice of the test case grid colored by the p initial condition from Equation ( |31| ). The inner box 
indicates the region where BASS’s solution was compared with EVA. 
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Figure 9. Slice of the test case grid colored by the pw initial condition from Equation ( |31 [ ) . The inner box 
indicates the region where BASS’s solution was compared with EVA. 


D. Spatial differencing scheme verification 

1. Linear analysis 


In this section, the results of verifying the implementation of BASS’s spatial difference schemes with the 
EVA tool will be presented. While not required to verify a numerical method, analysis of the accuracy of 
a scheme is often helpful in interpreting code verification results. The accuracy of both spatial differencing 
and time-marching schemes are commonly analyzed using Fourier analysis. For spatial differencing schemes, 
this approach consists of investigating the error associated with approximating the derivative of a single 
component of a Fourier series fk{x) = e ikx , where k is a real constant (the wavenumber) and i = \J — 1. 
Using the second-order central scheme 


D x (f(x,t )) = 


f(x + Ax, t) - f(x - Ax , t) 
2Ax 


to differentiate fk{x) gives 


D x (fk(x)) = 


f k (x + Ax) - f k (x- Ax) 


1 

2Ax 


2Ax 

( ik[x-\-A.x\ _i_ ik[x—Ax\ \ 

-he j 


= — [sin(fcAx)] e ikx , 


which is then compared to the exact value 


= ike ikx = * e ikx_ 

dx Ax 1 1 


The bracketed quantity in (34), 

( kAx )* = sin(fcAx), 

defined as the numerical wavenumber, appears when evaluating the relative error 

o,(A) - _ (tax)- 

e (kAx) — - kAx 

OX 


(33) 


(34) 


(35) 

(36) 


( 37 ) 
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associated with a spatial differencing scheme. Also of interest is p , the convergence rate of the relative error, 
which can be calculated by considering the slope of e(kAx) on a log-log plot, i.e., 


P = 


dlog(e) 
dlog(fcAx) 
dlog(e) dkAx 
dkAx dlog(fcAx) 

dlog(e) ; 
dkAx 
1 de 
e dkAx 


-kAx 


kAx. 


(38) 


The numerical wavenumber, relative error, and convergence rate of the spatial differencing schemes 
available in BASS are shown in Figure [lO] and the top and bottom of Figure [IT] respectively. The latter plot 
uses the number of points-per-wavelength for the abscissa. 



kAx 


Figure 10. Numerical wavenumber (kAx)* of the four spatial differencing schemes implemented in BASS. 


A few features of Figures 10 and 11 are worth mentioning. For the numerical wavenumber results 
in Figure [lOl the performance of each scheme can be judged by how closely it follows the curve corresponding 
to the exact value (kAx)* = kAx. Clearly, then, the explicit second-order and compact sixth-order schemes 
are the least and most accurate, respectively. The DRP and explicit sixth-order schemes make for a more 
interesting comparison, however: the DRP scheme outperforms the explicit sixth-order for higher kAx 
values, but it becomes difficult to distinguish the two schemes as kAx decreases. Inspection of the relative 
error plot in Figure 0 however, plainly shows that the DRP scheme’s error is actually higher than the 
explicit sixth-order for high points-per-wavelength (low kAx). The DRP scheme is an optimized form of 
the explicit sixth-order: it is essentially the explicit sixth-order scheme with modified coefficients that give 
better performance for moderate points-per-wavelength. 

The convergence rate curves clearly identify the order of the truncation error of each scheme: the explicit 
second-order scheme is second-order, the DRP scheme is fourth-order, and the explicit and compact sixth- 
order schemes are sixth-order. Considering the top and bottom plots of Figure 0 helps illustrate the 
difference between accuracy and order-of-accuracy. For example, despite its lower order-of-accuracy, the 
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Figure 11. Relative error and convergence rate of the four spatial differencing schemes implemented in BASS. 


DRP scheme is more accurate than the explicit sixth-order scheme for a considerable range of points-per- 
wavelength, and even outperforms the very accurate compact sixth-order scheme for a narrow band of ■ 
And though they possess identical orders-of-accuracy, the explicit and compact sixth-order schemes return 
widely different error over nearly the entire range of kAx shown in Figure El with the compact scheme 
beating the explicit by about an order-of-magnitude. 


2. BASS results 


As discussed above, the goal of code verification is to compare the convergence rate of the error of a PDE 
code to the formal order-of-accuracy of the code’s numerical method. For a spatial differencing scheme, this 
involves observing how the error behaves as the grid is refined. For the present test cases, each grid began 
as a uniform Cartesian grid extending from —0.5 to 0.5 in the three coordinate directions, and was then 
rotated, skewed, and subjected to sinusoidal perturbations. The coarsest grid in the series used 9 3 points — 
denser grids were constructed by uniformly refining the base grid through the addition of 1, 2, 3. . .points 
between each of the original grid’s, giving a grid series of 9 3 , 17 3 , 25 3 , 33 s , . . . , 201 3 points with approximate 
grid spacing 


Aa; = 


1 


(39) 


ranging from | to Aq. Figure 12 shows the 41 3 -point example from this grid series. 

For each BASS run, the RK7S time-marching scheme was used to march the initial flow to a final time 
level of t = 0.0005 with one time step. The choices of very small At and highly accurate time-marching 
scheme are meant to ensure that the error associated with the time-marching scheme is much lower than 
that of the spatial differencing. For the EVA calculation, a desired truncation error of 10“ 15 was used, which 
was easily satisfied, again thanks to the small time step size. 

The 1 2 e rror and the convergence rate of the I 2 norm, pi 2 , for the density flow variable, are shown 
in Figure 13 The quantity 2n/(kAx), an approximation of the number of points-per- wavelength, is used for 
the abscissa to facilitate comparison with Figure [TT] While not identical, Figures |TT] and [13] are remarkably 
similar. Like the linear analysis, the compact sixth-order and explicit second-order schemes are clearly the 


17 of Eg] 


American Institute of Aeronautics and Astronautics 



Figure 12. BASS verification test case grid with 41 3 points. 


most and least accurate, respectively. The DRP scheme outperforms the exphcitjdxth-order scheme for 
moderate points-per- wavelength, but then returns higher error for the high 2n/(kAx) values, as predicted 
by Figure |11| Many similarities are found in the convergence rate results, also: the compact sixth scheme 
initially shows higher-order behavior than the explicit sixth; the characteristic order spike is seen in the DRP 
curve, and the explicit second-order converges at its expected rate for a wide range of points-per-wavelength. 


Overall, it made little difference if the convergence rate was calculated using (27) or (30). Figure 14 
compares the two approaches for the l m ax norm of the pu error. The expected convergence rate is eventually 
observed for each scheme, and the curves are qualitatively similar. This indicates that the temporal compo- 
nent of the error is small relative to that of the spatial differencing scheme, and thus the assumption of ( 25 ) 
is a good one. 

Achieving asymptotic convergence is the most difficult for the compact sixth-order scheme. Figure [15] 
compares the compact and explicit sixth-order performance for all flow variables. Notice how significantly 
more points-per-wavelength are required for the compact scheme to achieve the expected p, a result not 
predicted by the linear analysis shown in Figure ED The cause may be related to the Giles non-reflecting 
boundary conditions used by the BASS code: any discrepancy between these BCs and the EVA solution 
will be propagated more readily by the compact scheme, as it requires a global matrix solve along each grid 
line. One would expect errors from the boundaries to become less significant as the number of grid points 
increases, which corresponds well with the behavior seen in Figure [To] 

All of the test cases discussed thus far have involved serial BASS calculations, i.e., the code used just one 
processor to compute the entire flow field. BASS, however, is a parallel code, and can decompose the solution 
domain into blocks that are distributed to and solved by many processors. To investigate the portions of 
the code responsible for the parallelization, the grid in Figure 12 was split into 8 3 = 512 blocks, and then 
run with the same initial condition, grid spacings, and time step sizes as the preceding test cases with 24 
processors. Figure [Tg] displays the grid’s block boundaries colored by block index number, showing how 
multiple blocks were placed inside the EVA domain to ensure any problems with BASS’s parallel routines 
would be noticed. 
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27t/ ( kAx ) 


Figure 13. p I 2 error and convergence rate for each of BASS’s spatial differencing schemes. pi 2 was calculated 
using (J 27 J . 




Figure 14. pu l m ax error and convergence rate for each of BASS’s spatial differencing schemes, comparing 
convergence rate calculation methods ( [27} and i pfij i. 
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Figure 15. I 2 error and convergence rate for BASS’s sixth-order schemes, with \27\ pi 2 calculation. 



Figure 16. Multiblock grid used for parallel test cases colored by block index. The extent of the EVA solution 
domain is outlined in black. 
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Figure |~iT| shows the density I 2 error and convergence rate results for each of BASS’s spatial differencing 
schemes (the equivalent of Figure [l3| . The smaller points-per- wavelength grids from the serial test case were 
not run because they did not meet the minimum points-per-block required by BASS. Like the previous test 
cases, each scheme’s convergence rate eventually attains its expected value. 




27r/(kAx) 


Figure 17. p I 2 error and convergence rate for each of BASS’s spatial differencing schemes. Parallel test case. 

The close agreement between the parallel and serial results in Figure [17] and Figure [l3| is unsurprising, 
since BASS’s numerical methods do not change when the code is run in parallel. The one exception is the 
compact sixth scheme, which uses an explicit sixth-order stencil at interior block boundaries. As described 
inp° this stencil has been tuned to give similar accuracy performance to the compact sixth scheme, and the 
results shown in Figure [Ts] confirm that this is the case. 


E. Time-marching scheme verification 

1. Linear analysis 

Here, the results of using EVA to verify BASS’s time-marching schemes are presented. As in the spatial differ- 
encing scheme verification section, Fourier analysis aids in the interpretation of the verification data, which, 
when applied to time-marching schemes, usually involves solving a simple ordinary differential equation like 


du 

— = —icon 
dt 


(40) 


with exact solution 


u(t) = uoe 


(41) 


where tt 0 = u(t = 0) and the time-marching scheme to be analyzed is used to solve (40). Following this 
procedure with time step size At for, say, Euler’s first-order explicit scheme 


u(x, to + At) 


u(x, t 0 ) + At 


dn 

dt 


x,t 0 


(42) 
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27 r/(kAx) 

Figure 18. pv I 2 error and convergence rate for the compact sixth-order scheme for the serial and parallel test 
cases. 


gives 


u(At) = uq + At (— icouo ) 
= (1 — icoAt) uo 
= G(ojAt)uo. 


(43) 


r< u(At.) 

u( 0) 


the ratio of the solution at the next and current time step, is called the amplification factor and 
is analogous to the ( kAx )* for spatial differencing schemes. Once a particular scheme’s G is known, the 
relative error associated with one time step (called the local error, Cf. 
expression 

^scheme (At) ^exact (At) 


36 p. 66]) can be found through the 


9(uAt) = 


^exact ( At ) 

= G(uiAt)/g(ujAt) — 1 


(44) 


where g(u>At) = e lu>At is the amplification factor for the exact solution (41). The convergence rate p for 
the error in (1 441) is found in a manner similar to a spatial differencing scheme’s, namely, 


1 d 6 . 

p = -- — —to At. 
9 dujAt 


(45) 


The magnitude of the amplification factor G for each of BASS’s schemes is shown in Figure [T9] and the 
corresponding relative local error and convergence rate curves for these schemes are shown in Figure |20| 
Similar to the (kAx)* plot of Figure [To} the accuracy of a scheme in Figure [19] is determined by how closely 
it matches the exact value of the amplification factor magnitude, which for (|41|) is 1. Figure 19 indicates 


the stability limit of the schemes, which is the largest value of to At where |G| < 1, i.e. , the amplitude of the 
solution does not grow after a time step. The amplification factor curves show, then, that the RK56 has both 
the highest accuracy and most restrictive stability limit of all of the schemes. The RK4L has approximately 
the same stability limit of the RK56, but dampens the solution considerably. The RK5L is clearly the least 
accurate scheme, but is much more stable than the RK4L and RK56. Finally, the RK67 and RK7S appear 
to possess the highest stability limits. 
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Figure 19. Magnitude of the amplification factor for each of the five time-marching schemes implemented in 
BASS. 




Figure 20. Magnitude of the relative local error and convergence rate for each of the five time-marching 
schemes implemented in BASS. 
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Figure [20] allows for a better comparison of the performance of the time-marching schemes. Like the 
corresponding spatial differencing error plot presented above, the quantity 27r/(wAf), the number of time 
steps per period of the solution in ( 41 ) is used for the abscissa. The error curves confirm that the RK56 
is a highly-accurate scheme, showing the sharp error “dip” characteristic of optimized numerical methods 
(e.g. the DRP spatial differencing scheme in Figure 11), and that the RK5L is by far the least accurate. 
The convergence rate curves show that the RK5L scheme eventually exhibits third-order behavior, while all 
of the other schemes converge at a fifth-order rate. The convergence rate of a Runge-Kutta’s local error is 
always one greater than its formal order-of-accuracy, indicating that the RK5L is second-order, and all other 
schemes are fourth-order. The RK4L scheme, however, is actually fourth-order for linear problems only. 
This scheme, the classic RK4 scheme (Cf. [36] p. 98, (235i)]) recast as a 2N Runge-Kutta method, drops to 
second-order when used to solve non-linear ODEs. 



Figure 21. Magnitude of the global amplification factor for each of the five time- marching schemes implemented 
in BASS after marching to ujT = 2n. 


Equation (44) shows how the error committed during one time step varies with At. Alternatively, one 
can calculate the error accumulated after marching to a constant final time level T = mnAt using n time 
steps for a ra-step scheme (the global error), which is 


$(n) 


^scheme ( ^~*) Inexact ( -^) 

Inexact (-^) 

[<?(*£)] 


ff(wT) 


- 1. 


(46) 


An expression for the global error’s convergence rate can be found by substituting uA t 
which eventually gives 

n d$ 

P $ dn' 


ujT/(mn) into (45 1, 


(47) 


Figures [21] and [22] show the magnitude of the global amplification factor, error, and convergence rate. Again, 
the RK56 and RK5L are seen to be the most and least accurate, respectively. The global convergence 
rate results show that each scheme converges at its formal order-of-accuracy. Also note how each scheme’s 
asymptotic range (the range of uAt the scheme converges at the expected rate) differs, with the RK56’s 
starting at the highest points-per-cycle value. Typical of low-order methods, the RK5L achieves its design 
order-of-accuracy more quickly than the other schemes. 
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Figure 22. Magnitude of the relative global error and convergence rate for each of the five time-marching 
schemes implemented in BASS after marching to c oT = 2n. 


2. BASS results 


To investigate the accuracy of BASS’s time-marching schemes, the procedure outlined in [B] was used with 
the initial condition found in (31 1 and Table [lj a single 201 3 -point grid with the same shape as the example 
shown in Figure |l2| and the explicit sixth-order spatial differencing scheme. BASS was marched to a constant 
final time level of t = 0.05 with 1 to 40 time steps, with the EVA tool provided the reference solution at the 
same final time level with a target truncation error of 10~ 12 , requiring a 26 th -order Taylor series. 

Similar to the preceding spatial verification results, an approximate frequency uj is used to plot the_error 
and convergence rates as functions of time steps per period of the disturbance. Here, the choice ui = ck was 
made, where c and k are the approximate maximum propagation speed and wavenumber from the spatial 
verification cases, respectively. 

Figure [23] shows the density I 2 error and convergence rate results for the temporal test case, with the 


convergence rate calculated from (27). The high-frequency portion of the plot (i.e., the left side) resembles 


the global error plot from the linear analysis, Figure [22] the RK56 is the most accurate, followed closely by 
the RK7S and then RK67, with the RK4L and RK5L showing significantly higher error. The results for the 
low- wavenumber runs, however, do not line up well with Figure |22[ with the exception of the RK5L, the 
schemes’ error curves flatten out, indicating that the error in the BASS calculation is no longer sensitive to 
the time step size. This, of course, is reflected in the convergence rate results. The RK7S and RK67 schemes 
only briefly attain their design convergence rates, and the RK56 never does. Also, the RK4L, a second-order 
scheme for non-linear problems, behaves like a fourth-order scheme before encountering the “error floor.” 


The density I 2 data from Figure 23 is reproduced in Figure ^4] but the convergence rate is calculated 


using ( 30 1 , which, unlike the spatial verification results, has a significant impact on the p data. The three 


fourth-order non-linear schemes now convincingly converge at a fourth-order rate, and the RK5L’s error still 
shows strong second-order behavior. The RK4L scheme’s error convergence, however, remains something 
between third and fourth-order. 


The RK4L’s convergence rate results in Figure 24 would seem to indicate that the present test case 
does not contain sufficient nonlinearity for the scheme to converge at its non-linear order-of-accuracy. To 
determine if this is the case, the RK4L scheme was run with an initial condition identical to that shown 
in (31) and Table [lj but with larger Gaussian perturbation amplitudes, specifically a = 0.1, u = 0.012, 
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Figure 23. p I 2 error and convergence rate for each of BASS’s time-marching schemes. Equation ( |27[ ) was used 
to calculate pi 2 . 




Figure 24. p I 2 error and convergence rate for each of BASS’s time-marching schemes. Equation ( |30| ) was used 
to calculate pi 2 . 
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27r/ ( uAt ) 


Figure 25. p I 2 error and convergence rate for BASS’s RK4L scheme with low- and high-amplitude Gaussian 
perturbations. Equation ( |30| > was used to calculate pi 2 . 


v = 0.08, w = 0.0075, p = jg-. The density I 2 results for this high-perturbation test case are compared 
to the original in Figure [25} As might be expected, the larger Gaussian amplitudes increase the I 2 error 
markedly, but also reduce the RK4L’s convergence rate to 2, the expected value for non-linear differential 
equations. 

The RK56’s Z max error exhibits an interesting behavior for the pw variable. Figure 26 compares this 
data with the corresponding 1-2 case. The pi 2 curve is fairly close to the expected fourth-order convergence, 
but p/ max is more like 3.5, a quality not seen in the other flow variables or time-marching schemes. The 


explanation may lie in the assumed form of the error used for calculating the convergence rate. For (30), the 
assumption 

e%,j ~ AAx f + BAt? (48) 

was made. An improved expression may be 


ejj ~ A At Ax? + BAt'j. 


(49) 


If (49 ) holds, one would expect the second term to dominate the error for most of the time steps in Figure [26] 
since the grid spacing is quite small. As the time step size is reduced, however, the two terms in ( |49[ ) may 
become closer in magnitude, which would tend to degrade the convergence rate calculation. This would 
especially be true for the RK56, the most accurate time-marching scheme, and the very sensitive Z max norm. 


V. Conclusions and Future Work 

In this work, the EVA approach to code verification was extended to verify the non-linear Navier-Stokes 
equations for the first time. Unlike the Method of Manufactured Solutions (MMS), the current state-of-the- 
art code verification technique, EVA does not require the addition of complicated source terms to a code’s 
governing equations, and thus the code itself. A survey of the relevant PDE code literature appears to show 
that the MMS source terms limit the technique’s popularity among researchers. Unlike the Method of Exact 
Solutions (MES), EVA maintains most of the flexibility in the form of the reference solution that MMS 
enjoys. 
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2ir/(uiAt) 


Figure 26. pw I 2 and £ max data for BASS’s RK56 scheme. Equation | |30| l was used to calculate both convergence 
rates. 


The EVA solution technique was explained briefly through its application to the linear advection equation, 
and then in some detail with regards to the Navier-Stokes equations. Pseudocode of the current EVA 
implementation were included and explained. The recurrence relations needed to construct the EVA solution 
are found in the appendix of this work. 

The EVA tool was then used to verify the viscous form of the NASA Glenn’s Broadband Aeroacoustic 
Stator Simulator (BASS), a high-order, parallel Computational Aeroacoustics code. After a Fourier analysis 
of BASS’s numerical schemes, verification results for the code’s spatial and temporal schemes were presented. 
For the most part, the results closely mirrored what was predicted by the Fourier analysis. Happily, the 
error from each of BASS’s schemes attained their design convergence rate, strongly indicating all are coded 
properly. 

The most significant limitation of the EVA method is the inability to specify boundary conditions along 
domain surfaces. In principle, it should be possible to use one of the spatial variables as the Cauchy marching 
direction instead of time, which would allow one to enforce a boundary condition at the cost of relinquishing 
control over the initial flow. This approach would require that the surface is not characteristic, and would 
not be applicable to solid wall boundary conditions. 

Another difficult aspect of code verification is determining the right set of grid spacing, time step size, 
and flow parameters that will allow the asymptotic range of the numerical scheme(s) under consideration 
to be found without extremely dense grids, many time steps, etc.. Assuming a more general form of the 
error than (28) for the order-of-accuracy calculation, and possibly refining grid spacing and time step size 
simultaneously, may be more flexible and robust. 
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Appendix 

The recurrence relation for the three-dimensional Navier-Stokes continuity equation used in this work is 


Qb-\-d-\- f -\-n-\- 1 g. 

dx b dy d dzf dt n+1 





j^b—a-\-d—c-\-f—e-\-n—m,y^ ^a+l+c+e+m ^ 

dx b ~ a dy d ~ c dzf~ e dt n ~ m dx a+l dy c dz e dt m 

Qb—a-\-d—c-\-f—e-\-n—rriy ^a+c+1 -\-e-\-m^j 

dx b ~ a dy d ~ c dzf~ e dt n ~ m dx a dy c+1 dz e dt m 

Qb—a-\-d—c-\-f—e-\-n—rriyj ^a+c+e+l+m^ 

dx b ~ a dy d ~ c dzf~ e dt n ~ m dx a dy c dz e+1 dt m 

Qb-a-\-d—c-\-f—e-\-n—m^j / 

dx b ~ a dy d ~ c dzf~ e dt n ~ m V 

Qa+l+c+e+m u 

dx a+1 dy c dz e dt m 

Qa+c+l+e+m v 

dx a dy c+1 dz e dt m 

ga+c+e+l+m w \ 

dx a dy c dz e+1 dt m ) 


( 16c repeated) 


and the corresponding recurrence relation for the x-momentum equation is 


gb+d+ /+n+l ^ 

dx b dy d dzf dt n+l 





f~jb—a-\-d—c-\-f—e-\-n—m ^ ^a+l+c+e+m,^ 

dx b ~ a dy d ~ c dzf~ e dt n ~ m dx a+1 dy c dz e dt m 

Qb—a-\-d—c-\-f—e-\-n—rriy ^a+c+l+e+m^ 

+ dx b ~ a dy d ~ c dzf~ e dt n ~ m dx a dy c+1 dz e dt m 

Qb-a-\-d—c-\-f— e-\-n—m yj ^a+c+e+l+m^ 

+ dx b ~ a dy d ~ c dzf~ e dt n - m dx a dy c dz e+1 dt m 

Qb—a+d—c+f—e+n—mo. ^a+l+c+e+m 

-I - 

dx b ~ a dy d ~ c dzf~ e dt n ~ m dx a+1 dy c dz e dt m 

Qb-a+d—c+f-e+n—m ^ ^ Qa+2+c+e+m u 

;i dx b ~ a dy d ~ c dzf~ e dt n ~ m \ 3 dx a + 2 dy c dz e dt m 

^ Qa+l+c+l+e+m v 

+ 3 8x a+1 dy c+1 dz e dt m 

1 Qa+l+c+e+l+m w 

+ 3 dx a+1 dy c d z e+1 dt m 

ga+c+2+m+m u ga+c+e+2+m u \ 

dx a dy c+2 dz m dt m dx a dy c dz e+2 dt m ) 


(50) 
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and the corresponding recurrence relation for the y-momentum equation is 


f +n+l y 

dx b dy d dz? dt n+1 


= -£ 



a — 0 

Qb—a+d—c+f—e+n—rriy^ 


^a+l+c+e+r; 


dx b ~ a dy d ~ c dzf~ e dt n ~ m dx a+1 dy c dz e dt m 

Qb-a-\-d—c-\-f—e-\-n—my Qa-\-c-\-l-\-e-\-m y 

dx b ~ a dy d ~ c dzf~ e dt n ~ m dx a dy c+1 dz e dt m 

Qb-a-\-d—c-\-f— e-\-n—m yj Qa-\-c-\-e-\-l-\-rriy 


dx b ~ a dy d ~ c dzf~ e dt n ~ 1 

Qb—a-\-d—c-\-f—e-\-n—m^j 


dx a dy c dz e+1 dt r 

Qa+c+l+e+m p 


dx b ~ a dy d ~ c dzf~ e dt n ~ m dx a dy c+1 dz e dt r 

Qb—a-\-d—c-\-f—e-\-n—m^j ^a+c+2+e+m^ 

v-. 


dx b ~ a dy d ~ c dzf~ e dt n ~ m 

^ ga+l+c+l+e+m, u 

3 dx a+1 dy c+1 dz e dt m 
1 d a+c+1+e+1+m w 


3 dx a dy c+ 2 dz e dt m 


3 dx a dy c+1 dz e+1 dt n 

Qa+2+c+m+m v 


+ 


d' 


a+c+e+2+m r . 


dx a+ 2 dy c dz m dt m dx a dy c dz e+2 dt m 


and, for the ^-momentum equation, 

Qb+d+f +n+l yj 

dx b dy d dzf dt n+l 


= -£ 



a— 0 

^b—a-\-d—c-\-f—e-\-n—rriy 


^a+l+c+e+m 


W 


dx b ~ a dy d ~ c dzf~ e dt n ~ m dx a+1 dy c dz e dt m 

Qb—a-\-d—c-\-f—e-\-n—nriy ^a+c+l+e+m^ 

dx b ~ a dy d ~ c dzf~ e dt n ~ m dx a dy c+1 dz e dt m 

Qb—a+d—c+f—e+n—m w Qa+c+e+l+m w 

dx b ~ a dy d ~ c dzf~ e dt n ~ m dx a dy c dz e+1 dt m 
dx b ~ a dy d ~ c dzf~ e dt n ~ m dx a dy c dz e+1 dt m 




Qb-a-\-d—c-\-f—e-\-n—m^j ^a+c+e+2+ra 


dt n 


dx b ~ a dy d ~ c dzf~ 

^ Qa+l+c+e+l+m u 

3 dx a+1 dy c dz e+1 dt m 

^ Qa+c+l+e+l+rriy 

3 dx a dy c+l dz e+1 dt m 

Qa+2+c+m+m w 

+ 


W 


3 dx a dy c dz e+ 2 dt n 


Qa+c+2+e+n 


W 


dx a+ 2 dy c dz m dt m dx a dy c+ 2 dz e dt 
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Finally, the energy equation recurrence relation is 


where 


Qb-\-d-\-f -\-n-\-lp 

dx b dy d dzf dt n+1 


= E 



a— 0 

Qb—a-\-d—c-\-f—e-\-n—my^ 


n 
m 

^a+l+c+e+m 


p 


dx b ~ a dy d ~ c dzf~ e dt n ~ m dx a+1 dy c dz e dt m 

£jb—a-\-d—c-\-f—e-\-n—rriy ^a+c+l+e+nip 

dx b ~ a dy d ~ c dzf~ e dt n ~ m dx a dy c+l dz e dt m 

Qb—a-\-d—c-\-f—e-\-n—rriyj ^a+c+e+l+m^ 

dx b ~ a dy d ~ c dzf~ e dt n ~ m dx a dy c dz e+1 dt m 

Qb— a-\-d— c-\- f — e-\-n—m p 

7 dx b ~° dy d ~ c dzf~ e dt n ~'‘ 

ga+l+c+e+m 

dx a+1 dy c dz e dt m 

Qa+c+l+e+rriy 

dx a dy c+1 dz e dt m 

Qa+c+e+l+m w 

dx a dy c dz e+1 dt r 

Qb+d+f+riQ 


- (7 - 1) 


+ 


Qb+d+f+nq, 


dx b dy d dzf dt n dx b dy d dzf dt 


Qb+d+n+nQ 

dx b dy d dz n dt n 


R ^ 

a—0 



Qb—a-\-2-\-d—c-\-f—e-\-n—m^j. Qa-\-c-\-e-\-rrip 

dx b ~ a+2 dy d ~ c dzf~ e dt n ~ m dx a dy c dz e dt m 

Qb-a-\-l-\-d—c-\-f—e-\-n—m^j ^a+l+c+e+m^ 


dx b ~ a+1 dy d ~ c dzf~ e dt n 

Qb-a-\-d—c-\-f—e-\-n—m^j 


71 dx a+1 dy c dz e dt ri 

ga+2+c+e+nip 


dx b ~ a dy d ~ c dzf~ e dt n ~ m dx a + 2 dy c dz e dt m 

Qb-a-\-d—c-\-2-\-f—e-\-n—m^j Qa-\-c-\-e-\-rrip 

dx b ~ a dy d ~ c + 2 dzf~ e dt n ~ m dx a dy c dz e dt m 

Qb— a+d— C+1+/— e+n — m ^ aa+c+l+e+m 

2 — 

dx b ~ a dy d ~ c+1 dzf~ e dt n ~ m dx a dy c+1 dz e dt m 

£jb—a-\-d—c-\-f—e-\-n—mg. ^a+c+2+e+mp 

dx b ~ a dy d ~ c dzf~ e dt n ~ m dx a dy c + 2 dz e dt m 


Qb—a-\-d—c-\-f—e-\-2-\-n—m^. 


9 


ia+c+e+m 


p 


+ 2 


dx b ~ a dy d ~ c dzf ~ e + 2 dt n ~ m dx a dy c dz e dt m 

Qb— a+d— c+f— e+l+n— ^a+c+e+l+mp 


dx b ~ a dy d ~ c dzf ~ e+1 dt n ~ m dx a dy c dz e+1 dt m 

Qb—a+d—c+f—e+n—m^j ^a+c+e+2+mp 


dx b ~ a dy d ~ c dzf~ e dt n ~ m dx a dy c dz e + 2 dt m 


( 53 ) 


( 54 ) 
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and 


Qb+d+f+nq, 

dx b dy d dzJ dt n 


b 

a— 0 




Qb—a+l+d—c+f—e+n—m u Qa+l+c+e+m u 

dx b ~ a+1 dy d ~ c dzf~ e dt n ~ m dx a+1 dy c dz e dt m 

f^b—a-\-d—c+l-\-f—e-\-n—my ^a+c+l+e+m^ 

dx b ~ a dy d ~ c+1 dzf~ e dt n ~ m dx a dy c+1 dz e dt m 

gb—a+d—c+f—e+l+n—m w ga+c+e+l+m w 

dx b ~ a dy d ~ c dzf ~ e+1 dt n ~ m dx a dy c dz e+1 dt m 

Qb—a-\-l-\-d—c-\-f—e-\-n—m^ ^a+c+l+e+m^ 

Q x b-a+ 1 dyd-c g z f-e Qfn-m Q x a g y c+l g z e gym 
Qb—a+l+d—c+f—e+n—m u ga+c+e+l+m w 

dx b ~ a+1 dy d ~ c dzf~ e dt n ~ m dx a dy c dz e+1 dt m 

gb—a+d—c+l+f—e+n—m^ ga+c+e+l+m w 

dx b ~ a dy d ~ c+1 dzf~ e dt n ~ m dx a dy c dz e+1 dt m 


Qb—a+l+d—c+f — e+n—m v Qa+l+c+e+m 

+ dx b ~ a+1 dy d ~ c dzf~ e dt n ~ m dx a+1 dy c dz e dt m 

Qb—a+l+d—c+f-e+n—m w ga+l+c+e+m w 

+ g x b-a+ 1 g y d-c g z f-e gfn-m g x a+ 1 g yC g z e g t m 
£jb— a+d— c+l+/— e+n— ^a+c+l+e+m^ 

+ dx b ~ a dy d ~ c+1 dzf~ e dt n ~ m dx a dy c+1 dz e dt m 

gb—a+d~c+l+f—e+n—m w ga+c+l+e+m w 

+ dx b ~ a dy d ~ c+1 dzf~ e dt n ~ m dx a dy c+1 dz e dt m 

Qb—a-\-d—c-\-f—e-\-l-\-n—my^ ^a+c+e+l+m^ 

+ dx b ~ a dy d ~ c dzf ~ e+1 dt n ~ m dx a dy c dz e+1 dt m 

Qb—a+d—c+f—e+l+n—rriy ^a+c+e+l+m^ 

+ dx b ~ a dy d ~ c dzf ~ e+1 dt n ~ m dx a dy c dz e+1 dt m 
+ 2 

Qb—a-\-l-\-d—c+f—e-\-n—rriy ^a+c+l+e+m^ 

dx b ~ a+1 dy d ~ c dzf~ e dt n ~ m dx a dy c+1 dz e dt m 

gb—a+d—c+l+f—e+n—m w ga+c+e+l+m ^ 

+ dx b ~ a dy d ~ c+1 dzf~ e dt n ~ m dx a dy c dz e+1 dt m 

gb—a+d—c+f—e+ 1+n—m ga+l+c+e+m w 

+ dx b ~ a dy d ~ c dzf ~ e+1 dt n ~ m dx a+1 dy c dz e dt m 


( 55 ) 
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